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^b ' We present a new numerical multiscale integrator for stiff and highly 

^~^ , oscillatory dynamical systems. The new algorithm can be seen as an im- 

fT^ ■ proved version of the seamless Heterogeneous Multiscale Method by E, 

Ren, and Vanden-Eijnden and the method FLAVORS by Tao, Owhadi, and 
Marsden. It approximates slowly changing quantities in the solution with 
higher accuracy than these other methods while maintaining the same com- 
Sr^ I putational complexity. To achieve higher accuracy, it uses variable meso- 

$H . scopic time steps which are determined by a special function satisfying 

moment and regularity conditions. Detailed analytical and numerical com- 
parison between the different methods are given. 



Abstract 



1 Introduction 

We consider numerical solutions of stiff and highly oscillatory dynamical sys- 
tems of the form 

^ = Mx), a;(0)=a;o (1) 



'Corresponding author. ylee@math . utexas . edu 
TengquistSmath.utexas . edu 



where the Jacobian of f^ has eigenvalues with large negative real parts or 
purely imaginary eigenvalues of large modulus. That is, the spectral radius 
is of the order. 

This imposes severe restrictions on the time steps. Resolving the e scale re- 
quires the time steps of a traditional direct numerical simulation (DNS) to be 
of order 0{e) or less. 

There are many numerical methods to approximate the solutions of ([l) with 
less computational complexity than 0{^) for 0(1) time intervals. Exponen- 
tial integrators or Gautschi type methods FlO HISI use an analytic form for the 
most significant part of the oscillatory solutions resulting in significantly less 
restriction on the time steps from stability and accuracy. Another method for 
highly oscillatory problems using asymptotic expansions in inverse powers of 
the oscillatory parameter ( [5] and references therein) has computational cost 
essentially independent of the oscillatory parameter. 

In this paper, we focus on the following two forms of the model problem 
which have scale separation. First, we consider the problem with explicitly 
identified slow and fast variables, 

§=/o(e,^), m^^o 

(2) 

— = , i]{0)=r]o, 0<e<Cl 

at e 

where i] is ergodic on some invariant manifold M{£,) for fixed ^. We also con- 
sider another problem 

§=/o(x) + ^, x(0)=xo, 0<6«1 (3) 

where the unperturbed equation 

dy fi{y) 



dt 



y(0) = 2/0 



is ergodic on some invariant manifold M (ya)- 

In (|2), 1] is called the fast variable because it has fast transient or highly os- 
cillatory behavior when the Jacobian of /i has negative real parts or all imag- 
inary parts. The slow variable £,{t) can be consistently approximated in any 
0(1) time by an averaged equation 



^=/(S):- / /o(S,77)dMS,r?) 



where /i(S, 77) is the invariant measure of 77 for fixed S. For more details, see 1191 
and II20I . In the case of l|3), it is often assumed for the analysis that there ex- 
ists a diffeomorphism from x to (^, 77) and this implies that there exist hidden 
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Figure 1: Schematics of HMM, MSHMM, FLAVORS and VSHMM 



slow variables in Q. The existence of slow variables for these problems mo- 
tivate the development of efficient numerical schemes for integrating the slow 
components of slow-fast systems without resolving all fast variables. 

In this paper, we focus on Heterogeneous Multiscale Methods (HMM) frame- 
work HJIZl that captures the effective behavior of the slow variables on the fly 
by solving the full scale problem in very short time intervals (see Figure [Tall. 
HMM does not require any a priori information about the effective force, and 
has a suitable filtering kernel to estimate the force by time averaging of the lo- 
cal microscopic solution. Because the effective force is independent of e and the 
fast variables, the time step for the slow variables which is called macro time 
step in HMM can be chosen independently of e. 

One of the variants of HMM, the seamless Heterogeneous Multiscale Method 
first introduced by Fatkullin, and Vanden-Eijnden in fS] and further devel- 
oped by E, Ren, and Vanden-Eijnden in [6J modifies HMM in that it does 
not require reinitialization of microscale simulation at each macro time step 
or each macro iteration step. In this strategy, the macro- and micro-models in 
(O use different time steps and exchange data at every step. The macroscale 
solver uses a mesoscale time step that is much finer than the one in HMM for 
the effective system, in order for the microscale system to relax and influence 
the macro scale (see Figure [Tbl l. We will here label the method MSHMM for 
mesoscale HMM in order to differentiate it from other methods called seam- 
less HMM (SHMM) for multi spatial and multi time scales without scale sepa- 
ration imEi. 

A similar technique can also be used to solve the system of the form (|3j 
without identification of the slow and fast variables beforehand. It was first 
noted by Vanden-Eijnden 11221 and later a variant was proposed and further 
developed by Tao, Owhadi and Marsden ||2T1 called 'Flow Averaging Integra- 



tors' (FLAVORS). It is based on the averaging of the instantaneous flow of the 
system with hidden slow and fast variables instead of capturing the effective 
force of the slow variables. By turning the stiff part on over a microscopic time 
step 5t and off during a mesoscopic time step h, FLAVORS obtains computa- 
tional efficiency (see Figure [Tel. 

In Section 121 we show that MSHMM and FLAVORS share a common char- 
acteristic in that they both approximate the effective behavior of Q and (|3j 
respectively by solving the problem with increased e values. 

The increase of the e value gives computational efficiency better than a di- 
rect approximation of the original problem but at the cost of reduced accuracy. 
The amplitude in highly oscillatory solutions, for example, is increased which 
is related to the increased e value. Because of this loss of accuracy, it is diffi- 
cult to generate higher order approximation of the effective behavior in both 
methods. 

The goal of the proposed method that applies to the more general formu- 
lation (|3]l is to increase accuracy by controlling the transient and the amplified 
oscillations while keeping the same computational efficiency and structure of 
implementation in the methods discussed above. And we call our proposed 
method a variable step size Heterogeneous Multiscale Method (VSHMM). To 
gain the control of the oscillations, the method uses variable mesoscopic time 
steps which are determined by a special function with a certain moment con- 
dition and regularity properties which is described in Section|3l Given a macro 
time step at which we want to sample the value of the averaged solution, the 
mesoscopic time step increases smoothly from fine one to a coarser mesoscale 
time step to obtain efficiency in computation. Once it reaches close to the next 
macro time step, the time step decreases again back to the original size and 
repeat this process for the next macro time step (see Figure ITdl ). 

VSHMM can also be used for many well-separated scale problems without 
using hierarchical iteration. Hierarchical iteration using the other multiscale 
methods - HMM, VSHMM and FLAVORS - has computational complexity 
which increases exponentially as the number of different scales increases HJ. 
Using the variable step size method, we can develop a new method whose 
complexity increases proportional to the number of different scales. Here, we 
focus on the two well separated scale problems and the new method for many 
scales will be reported in a forthcoming paper by the authors |,16J . The basic 
idea is to include different components of the force depending on the variable 
step size, from the full fc [x) for the shortest step size to only the slowest com- 
ponents for the longest step size. The intermediate step size will contain the 
intermediate to slow components of fe{x). 

As stated above, we have to mention that VSHMM requires scale separation 
and ergodicity of the fast variables. For stiff dissipative problems without scale 
separation, efficient methods exist such as implicit methods for small systems 
and Chebyshev methods for large systems. Here we study the application of 
VSHMM to dissipative problems with a potential application of VSHMM for 
concurrent multiscale problems in mind. 

This paper is organized in the following way. In Section |2l we review 



MSHMM and FLAVORS and show that they are equivalent in that they both 
solve a modified equation with increased e value. In Section |3l we propose a 
new method as an extension of MSHMM or FLAVORS to control the transient 
and the amplified errors and introduce higher order methods. In SectionlH we 
analyze the proposed method for dissipative and highly-oscillatory systems. 
In Section |5l numerical examples of dissipative and highly oscillatory systems 
are shown and also higher order method is verified. 

2 MSHMM and FLAVORS 

In this section, we review and compare MSHMM and FLAVORS. They are 
shown to share common characteristics except in modifications to the way 
time-stepping is implemented. 

The philosophy behind MSHMM lU is that we use different clocks for slow 
and fast variables. It requires identification of slow and fast variables in ad- 
vance and is applicable to l|2]|. If we denote the micro and mesoscopic time 
steps by 5t and h, an explicit first order MSHMM solves for rj first, 

^"+' = r;" + -/i(r,77") = 11- + ^fiicvn 

e e' 

with e' = e4^ and it uses the information from this calculation for the evolution 

dr 

ofC, 

This is a consistent approximation to the model problem with e modified to the 
increased value e' = ct^ > e. 

OT 

In 11211 , Tao et al. propose another method based on the averaging of the 
instantaneous flow of the system, which is called FLAVORS. It turns on and off 
the stiff parts to capture the effective flow of the slow variables. It solves the 
full problem (|3j with the stiff part with a micro step St, 



x" + Jf(/o(x") + 
It then uses a mesoscopic time step h without the stiff part 

If we consider this as one single step, it is 

a:"+^ = x-+Uh{x-) + /i/o (x- + 5t Uix") + ilip.\\\+i^^±J^f,{j:-) 

/i(2:") 



x" + iSt + h)[foix-+**) + 



with an e value increased to e' = ^^j^e and x"~^** such that 
/o(x"+**) = j^ {Stfoixn + hfoix^+*)) . 

Hence, it solves the model problem with modification to e and a minor differ- 
ence in the /o term. 

If we apply FLAVORS to (|2), it becomes much clearer that MSHMM and 
FLAVORS share common characteristics. As before, we use the explicit first 
order Euler method for each step for ^ and i]. We have 

and 

The evolution of ^ can be represented in a compact form, 

=r + (st+h) {ofoicvn + (1 - e)/o(r+*,^"+')) (4) 

=r + (St + h) {0fo{c\vn + (1 - o)fo{c\r^')) + o{6th) 

where 0=^. 
If we now choose 

St in MSHMM = Stm FLAVORS 

and 



h in MSHMM = 5t + h in FLAVORS, 



then they both are first order approximations (with slightly different fact that 
FLAVORS uses the 9 method for the slow variable) to the following modified 

'■ ■ h_ 

St' 



equation, with increased e' = ^Tj^e = (1 + 0^)^: where a ~ — 



^J=m,v), o<t<T 

1 1 - (5) 

For simplicity, the reduction in the overall processes of time integration, 
which is a, is called the savings factor. The savings factor gives information 
for computational efficiency of MSHMM or FLAVORS but it is not the actual 
computational efficiency of the methods. For FLAVORS, as an example, if the 
same order integrators are used for the forces with and without the stiff part, 
the number of function evaluations of the force terms for direct numerical sim- 
ulations (DNS) and FLAVORS using the same micro time step 5t are given by 



[1 + a] and 2 for the time 5t + h where [•] is the ceiling or the smallest integer 
function. Therefore, the computational efficiency of FLAVORS over DNS is 

For better computational efficiency, it is obvious to use larger a values. But 
arbitrarily large a values do not guarantee the convergence of the methods to 
slow variables. In f^T), the relation between 5t, h and e is analyzed for conver- 
gence: 

^t^ , c St 

-— -^h + St-^ —. (7) 

e^ e 

Using a instead of h, we can check 

(a + l)e < 1 (8) 

for the convergence of FLAVORS. 

To see the effect of the increased e' = {l+a)e value, we compare the solution 
^{t) of I© with the solution ^{t) of © for 0(1) time t. Let S(t) be the effective 
solution of l|2|. Because lO has the same invariant measure of l|2ll, S(f) is also 
the effective solution of l|5|. For the averaging error, it is normally expected to 
have the following type of error bound (see Iiln7nl9il for example) 

\\m~^moo<ce'^ (9) 

for constants a > and C which is dependent on time t and independent of e. 
Similarly, the averaging error of ^ by S is given by 

iie>)-s(t)ii^<c(i + «)"£'' 

which implies the following amplified error due to the increased e' value. 

\\m-m\.o<c{i+are\ 

If we want more computational efficiency, the savings factor a = ^ should 
be larger. However then we lose accuracy because of the amplified oscilla- 
tions. The key feature of the proposed method is to decrease O ((ae)°) term to 
0(e") which is independent of the savings factor, a, while keeping the same 
computational efficiency depending on a. Because we are looking for effective 
behavior of systems with 0{e) perturbations, the diminished oscillation and 
fluctuation to 0{e) is accurate enough to approximate the effective behavior of 
the slow variables. 



3 A Variable Step Size Mesoscale HMM (VSHMM) 

In this section, we propose a variable step size Heterogeneous Multiscale Method 
(VSHMM) which controls the transient and the amplified oscillations of MSHMM 



and FLAVORS while maintaining the computational complexity and general 
structure of the methods. The new method is a modification of MSHMM and 
FLAVORS. The key feature of the new method is to use variable sizes of meso- 
scopic time step. 

In dissipative problems, when all eigenvalues of the Jacobian of /i or the 
real part of them are negative, transient behavior of the fast variable to the 
quasi stationary state contributes a significant part of the error [9]. If e is mod- 
ified to a greater value, then the fast variables change slower, resulting in an 
error that remains in the quasi stationary solution after the transient. There- 
fore, it is necessary to use small e values at the beginning of each macroscopic 
time step to guarantee that the fast variables relaxed rapidly. 

In highly oscillatory problems, when the eigenvalues of the Jacobian of /i 
are imaginary, as we mentioned in the previous section, the increased e ampli- 
fies the oscillations and this dominates the error which can be controlled using 
higher order methods. If there is no a priori identification of slow variables, 
the only possible way to control this amplified error is to use a smaller e value 
locally which requires a finer time step. 

Our method, VSHMM, reconcile these contradictory situations by introduc- 
ing time dependent e values. At the beginning and the end of each macro time 
step, it uses very fine mesoscopic time steps to overcome problems such as the 
delayed relaxation of the fast variables in dissipative problems and amplified 
oscillations in highly oscillatory problems, while using coarse mesoscopic time 
steps at all other times to save computational complexity (see Figure[l]for com- 
parison of time stepping with the other method). Once the system has evolved 
to the next macro time step, we iterate the same process. 

Therefore, by using the variable mesoscopic step sizes, we expect to ob- 
tain a more accurate approximation of the slow variables than MSHMM and 
FLAVORS, only after macro time steps. We emphasize that there is no explicit 
macro time stepping in the new method but we use the calculated values only 
at the specified macro time intervals because other values are not guaranteed 
to give less amplified errors. In the highly oscillatory problems, if we want to 
have the same savings factor as MSHMM and FLAVORS, the intermediate val- 
ues of the new method between two macroscopic sampling time becomes more 
oscillatory than MSHMM or FLAVORS (see Section|5]for numerical examples). 
This is because e is modified for a larger value than the modified e of MSHMM 
or FLAVORS to compensate the loss of efficiency at the beginning and the end 
of each macro time step. 

3.1 Description of the new method 

In the proposed method, the mesoscopic time step h for MSHMM or FLAVORS 
is described as a time dependent function. Let K E C|((0, 1)) with a compact 



support in (0, 1) such that 



/ K{t)dt = 1, 
Jo 



(10) 

/o 

d'-A'(i) 



dr 



= 0, r = 0,1, ...,g fori = 0,1. (11) 



For a given macro time step AT, the time dependent mesoscopic time step h{t) 
is given by 

h{t) = adtKAT,qi^AT.qit mod AT)), (12) 

for a savings factor a > 1 when Kat.q is a rescaled version of K, 



AT ^AT' 



KAT,q ^l^Kil^) 



and QAT,q{t) is the antiderivative of Kat.q with 6at,<j(0) = 0. 

It can be easily verified that A'at,? satisfies the following moment and reg- 
ularity conditions 



KAT.qdt = AT, (13) 

d'-KAT,git) 



dt^ 



0, r = 0,l,...,q, t = 0,Ar. (14) 



Because can be seen a special case of Q, we describe the proposed 
method for the case of l|3j. 

ALGORITHM - one macro time step integration of VSHMM 

Let x" be the numerical solution to l|3]l at t = i" := n/ST with a savings factor 
a. 

1. Integrate the full system for 5t to resolve the fast time scale 

x{e + 5t) = $^,i(i") 

where $^j is an integrator of dx/dt — fo{x) + fi{x)/e for St. 

2. Update time 

t = e + St. 

3. Integrate the system without the stiff part with mesoscopic time step h{t) 

where $°(f) is an integrator of dx/dt = fo{x) for h{t). 

4. Update time 

t = t + h{t). 



5. If time is at the macroscopic time points, sample the solution 

6. Repeat from 1 for the next macro time step integration. 

The mesoscopic step size h{t) is very small when the simulation starts and 
smoothly increases. Once it reaches ■^, it decreases smoothly back to the small 
value again as t — > AT. Equivalently, the new method solves the system with 
small e value and the e value increases smoothly to accelerate the computation 
and returns to the original small value for the next coarse step. 

In VSHMM, the ratio between the mesoscopic and microscopic time steps 
is not constant. The following proposition illustrates that for a given savings 
factor a in lO, the computational efficiency of VSHMM is same as the case 
when the mesoscopic time step is constant with the same savings factor. 



Proposition 1. 



AT L 6t "'■ 



Proof. This is a simple application of change of variables. Let s = QAT,q{t), 
then 



\tJo St AT Jo ^^-'^^ 



ds 

a 



AT Jo St AT Jo •'"'KAT.qit) 

from the moment condition of iiTAT,?- CH 

3.2 Higher Order Methods 

If there is a priori identification of the slow and fast variables, MSHMM may 
be implemented with higher order methods and this can also be applied to 
VSHMM. Without identification of slow and fast variables, it is not easy to 
implement a higher order scheme for the slow variables. We show that for 
VSHMM a second order mesoscopic integrator, for example, an explicit Runge- 
Kutta method, gives quadratic decrease of errors with an additional error term 
which can be ignored in comparison with the dominating error. 

Here we describe a second order approximation. We integrate the full sys- 
tem with higher order numerical method for St to resolve the fast time scale, 

x(0=$|,(i°). 

Then we use the second order explicit Runge-Kutta method for /o() part with- 
out /i part in (|3j, 

r = S:{t)+'^fo{m) 
x{t + h{t)) = x{t) + h{t)fo{x*) 

10 



In MSHMM and FLAVORS, the amplified error dominates other error terms 
from mesoscopic and microscopic integrators. Because MSHMM and FLA- 
VORS cannot control this amplified error, it is difficult to see the effect of the 
higher order mesoscopic integrators. With VSHMM which controls the ampli- 
fied error, higher order mesoscopic integrators can be verified (see Figure |4] in 
Section|5]for a numerical result). 

4 Analysis 

We analyze VSHMM for highly oscillatory problems. First, we start with a re- 
view of the dissipative case and address the importance of the rapid relaxation 
of the fast variables at the beginning of the simulation. The following result for 
the dissipative problem is from flSll . 

Theorem 4.1. ilTgl f Assume that for fixed ^ of ([2J, 77 has a unique exponentially 
attracting fixed point, uniformly in ^. Specifically we assume that there exists p and 
a > such that, for all S, and all 771, 772, 

/i(C,p(O)=0, 

{fi{£.,Vi)~ fi{£.,V2),Vi -m) < -a\rii -772p. 
Also assume that there exists a constant C > such that 

\k{Ul)\<C, |V,/o(^,77)|<C 

|V,/o(e,^)|<C, \ii{x)\<C 



and 

If^{t) is the solution to 



|V?7(a;)| < C. 



jE = f{E,p{E)), E{0)=m, 

Then there are constants M, c > such that 

|^(i) - S(i)P < ce*"(e|77(0) ~ p(C(0))p + a 

This theorem indicates that the 0(1) error in ?/ for the first time step may 
give an ©(e) global error. If the e value does not change at the beginning of sim- 
ulation to guarantee that 77 relaxes close enough to p(C(0)) and then increases 
to a larger value after relaxation, for example ©(^/e), to get computational effi- 
ciency, the global error is still ©(e). Therefore, it is important to have a small e 
value at the beginning of simulation at each macro time step in VSHMM (see 
Figure|2]in Section|5]for a numerical example). 
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4,1 Highly Oscillatory Case 

In the highly oscillatory case, the effect of variable mesoscopic time integration 
is analyzed and we show that the error is of order e independent of a, which is 
significantly less than MSHMM and VLAFORS for large a values. 

Instead of regarding the new method as solving with an increased e, we can 
rescale time resulting in multiplication by a factor in the equation for the slow 
variables. Using the same procedure in Section |2l it can be verified that the 
explicit Euler version of the proposed method for ^ for < t < T with ([T2] l is 
equivalent to solving the following modified problem, 

7 rji 

-e =(1 + aKAT,M^ + a)t))fo{^, ri), 0<t< 



dt ^ ' 1 + a (^5) 

dt e 

where Kj^T.qit) satisfies | [T3] | and (|T4). Note that in the formulation above, we 
do not have Q^^{t) as an argument of K^T,q while the mesoscopic step sizes 
are given by lfT2|l which is 

hit) - a/vAT,,(e-i(t)). 

In many oscillatory situations, /^ assumes special forms such as fe{t) = 
fe{t,t/e) which are periodic in the second variable. We hypothesize that the 
effective force of the system can be defined by 



fit) = lim 

(5->0 



1 f'+^ 

lim^ / fe{T)dT 



as in ill and 0. 

Based on this hypothesis, the averaging error (|9) has a = 1 in highly os- 
cillatory problems which is ©(ae). The next theorem shows the effect of the 
rescaled system using the time dependent mesoscopic rescaling function. The 
rescaled system approximates the averaged system with an C(e) error term 
which is independent of e. 

Theorem 4.2. Let (^, rj) and (^, fy) be the solutions to and l fI5l > respectively with a 
savings factor a and K^r.q satisfying (I73l > and (Il4l > and the fast variables are periodic. 
Further assume that V xfoix, y) is hounded independently ofe. Then 

r— 1 

where Ci,i — 1,2,3, are constants independent ofe and a. 

Proof. We use notation /C(i) to denote (1 + aA'AT,g((l + a)i)) to simplify the 
argument. First, for fixed ^ and ^, i] and fj have the same invariant measure. If 
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we denote the averaged solutions to l|2ll and ((TSl l by S and S respectively, they 
satisfy 

" " "' < t < AT (16) 



and 



/C(t) satisfies 



AT 
" l + o 



d^^-^"(^)' 



"^ S = /C(f)/(S), 0<t< ^^ 



dt 



AT 

" l + o 



1 + a 



(17) 



]C{t)dt = / (1 + a/\AT,q((l + a)t)) dt 
lo Jo 



AT 



1 + a 1 + a Jq 

=AT 



I K^T,q{s)ds with s = (1 + a)t 
Jo 



(18) 



For flTt , by using dlSl and change of time r such that dr/dt — K[t), we can 
verify that 

~ AT 

1 + a 
Now we compare S(t^) and ^(^^)- First, 



e 



AT 
l + a 



AT 

1 + c 



^W/o(l,'7)di- 



Letgo(^,Vf) •= .foi(,it):V{t)) and /(t) = J go{t, s)ds where go is 1 -periodic in 
the second variable. Then 



^ f AT 



1 + a 



AT 

■ l + o 



K,{t)f{r)dt. 



Therefore we have 

AT 



e 



1 + a 



^ f AT 



1 + Q' 



AT 

■ 1 + c 



IC{t)gi{t,t/e)dt 



where 



gi{t,t/e) = go{t,t/e)^f{t). 



Using Lemma fOl we prove the theorem for the case when .go(^, t/^) is periodic 
in the second variable. D 



Lemma 4.3. 



AT 

1 + c 



ICit)gi{t,t/e)dt 



< Cie + C2Y, 



— a e ^ 



r=l 



AT 



ir-1 



C: 



(ae) 



9+1 



AT« 



w/zere Ci, i = 1, 2, 3 are constants independent of e, a and AT. 
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Proof. Partition the interval, (0, ^^), into N uniform subintervals, (i„, i„+i), n = 
0,1,...,N- 1 such that 

AT 

to ~ 0, tjv 



1 + a 



and 



\tn+l — tn\ — e. 



This condition requires that N = (i+Z)e ■ 

For t„ <t< tn+i,n = 0,1, ...,N— 1, using Taylor series expansion of gi(-, •) 
in the first variable at i = i,i+i/2 = -^^^-| — -, we have 

gi{t,t/e) = 9i{tn+i/2,t/e) + digi{tn+i/2,t/e){t - f„+i/2) + O(e^). 

||(?i5l|oo is bounded and independent of a and e. The second term gives e 
order term and after integration on [t„,t„+i], it becomes e'^ order. There are 
N = ,-1^-'", intervals, therefore, we have 



AT 

l + a 



N-1 



IC{t)gi{t,t/e)dt^ J2 



n=a 



^{t)9i{tn+i/2-,t/e)dt 



-0{ATe) 



We further analyze the first term on the right hand side using integration 
by parts, which gives 



N-l 



n=0 



tr.+i r'^+ 



^ e/C(t).gW(i„+i/2,iA) -/ e/C'(i),9W(t„+i/2,tA)di (19) 



where ^[^l (t„+i/2, s) isanantiderivativeof gi(t„+i/2, s) such that /(7[^l(t„_|.i/2, s)ds 
0. From this mean zero condition ^I^l (in+i/2i s) is periodic in s. 
For the first term of ((T9ll , after rearrangement of the summation, 

n=0 t„ I e 

- ^ /C(t„) (5['l(Wl/2, f) - 5W(in-l/2, ^)) 
n— 1 ^ ''^ 



-/C(io)gW(ii/2,-) 



(20) 



Using Lemma|13]with the facts that ||A:||oo = 0{a) and A^ = 0{^), we 
estimate the second term of l|20|l . 



N-l 



e J2 ntn) ( .9^(^1/2, -) - 5W(i«-i/2, -) ) = 0(AT6). 
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For t = to and t^, K{i) = 1, and we have 0{e) for the first and last terms of 
®. 

For the second term of l fT9] l, we do the same procedure except that now we 
have 

!l/c'(t)lloo = o(|^) 

/C'(t)=0 fort = io,iw, 

and 

.9W(in+i/2,s)-.gW(in-i/2,.s) = 0(e). 

After integration by parts, we have 

n=0 ^ *" ^ n=0 •^*" 

where g[^l(t„+i, s) isanantiderivativeof .g[^l(i„+i,s) such that /5[^l(i„+i/2, s)ds ~ 
0. Now Loo estimate of the right hand side is 

3 2 2 

^ AT^^ ^ AT ' 

If d^lC{t)/dV = for t = h,tN, and r < q, we can repeat the same procedure 
and this proves the lemma. D 

Lemma 4.4. 

5™ (Wi/2, s) - ff''' (in-i/2, s) = 0(e) 

Proof. For t„_i/2 < t* < i,i+i/2. 



where constant C(i*) is given by 



(i,,s)= / gi(i*,T)dr - C(i,) 





C(**) = / / gi{U,T)dTds. 
Jo Jo 

Now 

5'^'(iri+l/2,s)-.9'^'(in-l/2,s) == / (ffl (t„+l/2 , t) - .91 (t„-l/2, t)) dr 

"'0 

+ C'(^ri+l/2) - C{tn-l/2) 

Using 

9l{tn+l/2,r) = 9l{tn-l/2,T) +0(e), 

we have 

(.9i(Wi/2,t) -ffi(^n-i/2,'r))rf'r ^0{e) 



15 



Also for the constant terms, 

C(in+i/2) - C'(^n-i/2) = / / {9iitn+i/2, t) - gi{tn-i/2, t)) drds 

Jo Jo 

=0{e) 
This proves the Lemma. D 

5 Numerical Experiments 

In this section, we apply VSHMM to dissipative and highly oscillatory prob- 
lems. The result of VSHMM is compared with MSHMM and FLAVORS. The 
convergence of the second order scheme is also verified for (|3]l which does not 
have a priori identification of slow and fast variables. As the last test problem, 
VSHMM is applied to a stellar orbit problem |3 13. 14]. The comparison of 
three methods, MSHMM, FLAVORS, and VSHMM shows that VSHMM cap- 
tures slow variables with higher accuracy than MSHMM and FLAVORS. 

In choosing K for the variable mesoscopic time step, the ioo norm of K 
plays more important role than the regularity conditions at the boundary. The 
error depending on the regularity of K has the form of powers of (ae) and the 
restriction ||8ll for choosing a shows that iiq > 1, then this error is small. On the 
other hand, if || A'||oo is too large, then instantaneous value of h{t) + 6t at some 
time violates the condition ^ and this violation deteriorate the convergence 
of VSHMM. In all numerical examples, we use K{t) = (1 + cos(27r(t - 1/2))) 
which has 1 1 A' 1 1 oo = 2 and q = I. Also, the micro steps for all methods are same. 
Therefore, the computational efficiency between DNS, MSHMM, FLAVORS, 
and VSHMM follows ©. 

5,1 Stiff Dissipative Case 

We begin with a dissipative example to show the effect of resolving the tran- 
sient behavior accurately. 

The following problem has explicit form of slow and fast variables, 

d £,-1] 
dt'^=^- 

with initial value (2:(0), 2/(0)) = (-1,1). 

We can verify that the averaged equation is given by 

f = i+s, m = -i (22) 

Figure |2] shows the numerical result from the various methods for e = 
2 X 10"'' and a = 100. For VSHMM, MSHMM and FLAVORS, a fourth or- 
der Runge-Kutta method is used for full scale integrator while a second order 



(21) 
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-DNS 
MSHMM/FLAVORS 
VSHMM 




Figure 2: Dissipative case l(2T|l . Plot of the slow variable £, of l|2T)l computed 
by DNS, MSHMM, FLAVORS and VSHMM. e = 2 x IQ--* and a = 100 for 
VSHMM, MSHMM and FLAVORS. 



Runge-Kutta method is used for mesoscopic integration. In this case, MSHMM 
and FLAVORS produce errors greater than the error from VSHMM because of 
the error from inappropriately resolved transient behavior at the beginning. 
VSHMM approximates the slow variables more accurately than MSHMM and 
VSHMM at coarse time points with a time step AT = 0.2. 

5.2 Highly Oscillatory Case 

The next three ntmierical experiments are highly oscillatory problems where 
the eigenvalues of the Jacobian of /i are all imaginary. The first two problems 
are expanding spiral problems. In the first case, the fast variable has a fixed 
angular period while in the second problem, the fast variable has variable pe- 
riods depending on the slow variable and fully nonlinear. The last problem is a 
well studied system taken from the theory of stellar orbits in a galaxy l3l lT3lll4l . 

5.2.1 Constant Angular Period Case 

The first oscillatory case is an expanding spiral problem for complex x with 
constant angular period. The equation is given by 

dx X 5Real(x)x ix ^ , . 

— = - + rT^ + -^ 2:eC (23) 

at 4 \x\ e 

x{0) = 1. 
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(a) global 
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Figure 3: Constant angular period case l|23l l. Plot of the slow variable ^ = |x| of 
igi computed by DNS, FLAVORS and VSHMM. Global (left) and local (right) 
representations of each solution in [0,3]. e = 3^ and a = 50 for VSHMM and 
FLAVORS. 



The slow and fast variables are ^ 
be verified that 



and 1] ~ arg(a;) respectively and it can 



e 



dt 4 
di] _1 

'dt ~1 
(e(0),rKO))=(l,0) 



5^ cos(r?) 



(24) 



which has constant periodic force for £_. 

Figure [3] shows the slow variable over time with a locally magnified plot 
in the neighborhood of f = 4. As shown in Figure 13b [ the proposed method, 
like FLAVORS, has amplified oscillations except in the neighborhood of the 
specified macro time interval. Once it reaches the neighborhood of the macro 
time interval, the new method has less oscillations and converges faster to the 
averaged solution. 

Next we verify the order of accuracy for the second order mesoscopic inte- 
gration methods. In VSHMM, there are several error terms from various fac- 
tors - order of accuracy of each integrator and regularity conditions and the 
Loo norm of K for variable stepping. To see the second order behavior, other 
error terms must be well controlled to be much smaller than the error of the 
second order integration. For this purpose, we choose a = 52.67 which was 
obtained from numerical tests. Figure H] shows the errors of the first and sec- 
ond order methods with the analytic effective solution. For sufficiently small 
average mesoscopic step size, aSt, the first order method shows linear decrease 
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Figure 4: Errors of the first and second order VSHMM for | [251 l. e = 10 "^ and 

a == 52.67. 



of error as expected while the second order method shows quadratic decrease 
of error for relatively large mesoscopic step sizes. The error of the second order 
method becomes flat for much smaller mesoscopic step sizes because the error 
is dominated by the averaging error which is of order e. 

5.2.2 Fully Nonlinear Case 

The second test problem is also an expanding spiral problem in C^ where two 
spirals are coupled. This system is more general than the previous oscillatory 
problem in that it is not angular periodic with nonlinear /i . For a fixed £_, the 
periodicity of each component of ?/ depends on ^ and the oscillation of ^ comes 
from rji and 772 simultaneously which has an irrational initial ratio. 



(25) 



The slow and fast variables are given by ^ = (|a:|, \y\) and?? = {axgi{x) , a.rgi{y)) 
respectively. Figure |5] shows the result of VSHMM and FLAVORS for ( |25t with 
e = 5 X 10^'' and q = 50. A fourth order Runge-Kutta method is used for micro 
step simulation and a second order Runge-Kutta method for mesoscopic time 
step simulation. VSHMM captures the correct slow variable on the macro time 
interval, AT = 0.6. In Figure|3j it is clear that averaging the FLAVORS solution 
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(a) global 



(b) local 



Figure 5: Fully nonlinear case l(25]| . Plot of the first slow variable ^i = |a;| 
computed by DNS, FLAVORS and VSHMM. e = 5 x lO^^ and a = 50 for 
VSHMM and FLAVORS. 



improves the approximation of the effective solution. Figure |5] shows that this 
is not always the case. 



5.2.3 Stellar Orbit Problem with Resonance 

The last numerical experiment is a well studied system taken from the theory 
of stellar orbits in a galaxy 



r^ + a ri =£r2, 
r2 + fe^r2 =2erir2 

where ri is the radial displacement of the orbit of a star from a reference cir- 
cular orbit, and r2 is the deviation of the orbit from the galactic plane. Here 
t is actually the angle of the planets in a reference coordinate system. Using 
an appropriate change of variables, the system can be written in the following 
form O 





1 
e 


/ a 0\ 

-a 

6 

\0 -b o) 


x + 


/ \ 

xj/a 



\2xiX2/b) 



x(0) 



/1\ 





X e 



(26) 



with a ~ 2, and 6 = 1. In |2H3l, it is verified that in the case of a = ±26, the 
system is in resonance and has three hidden slow variables ^^ : R"' ^- M, i = 
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Figure 6: Stellar orbit problem l|26l l. Plot of the three slow variables (|27} com- 
puted by DNS (real line) and VSHMM (marked with cross, square and circle), 
e = 10-4 and a = 100. 



1,2,3, are given by 



6 
6 



(27) 



^xixl 



2x2X32-4 



xix\ 



The resonance of oscillatory modes generates lower order effects, that are cap- 
tured by VSHMM. In Figured we present a numerical result of our method for 
e = 10-4 and a = 100. 



6 Conclusions 

We have presented a new multiscale integrator VSHMM for stiff and highly 
oscillatory dynamical systems. It controls the transient and the amplified os- 
cillations of MSHMM and FLAVORS while preserving the computational com- 
plexity and general structure of these methods. This results in an overall higher 
accuracy. The main idea of the error control is to use variable mesoscopic step 
sizes determined by special functions satisfying moment and regularity con- 
ditions. The proposed method is restricted to ordinary differential equations 
with two scales. Applications to stochastic differential equations and problems 
with more than two scales will be reported in a forthcoming paper by the au- 
thors ini. 
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